Modelos de Nichos Ecológicos

Carga de paquetes

# Colección de paquetes de Tidyverse
library(tidyverse)

# Estilos para ggplot2
library(ggthemes)

# Paletas de colores de RColorBrewer
library(RColorBrewer)

# Paletas de colores de viridis
library(viridisLite)

# Gráficos interactivos
library(plotly)

# Manejo de datos vectoriales
library(sf)

# Manejo de datos raster
library(terra)

# Manejo de datos raster
library(raster)

# Mapas interactivos
library(leaflet)

# Acceso a datos en GBIF
library(rgbif)

# Datos geoespaciales
library(geodata)

# Modelado de distribución de especies
library(dismo)

library(rJava)

Obtención de datos de presencia

# Nombre de la especie
especie <- "Bradypus variegatus"

# Consulta a GBIF
respuesta <- occ_search(
  scientificName = especie, 
  hasCoordinate = TRUE,
  hasGeospatialIssue = FALSE,
  limit = 10000
)

# Extraer datos de presencia
presencia <- respuesta$data

# Guardar los datos de presencia en un archivo CSV
write_csv(presencia, 'presencia.csv')

# Leer los datos de presencia de un archivo CSV
presencia <- read_csv('presencia.csv')

Warning: One or more parsing issues, call `problems()` on your data frame for details,
e.g.:
  dat <- vroom(...)
  problems(dat)

Rows: 6619 Columns: 173
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr  (122): scientificName, issues, datasetKey, publishingOrgKey, installati...
dbl   (28): key, decimalLatitude, decimalLongitude, crawlId, taxonKey, kingd...
lgl   (17): isSequenced, isInCluster, acceptedNameUsage, locationRemarks, ta...
dttm   (5): lastCrawled, lastParsed, dateIdentified, modified, lastInterpreted
date   (1): georeferencedDate

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

# Convertir a conjunto de datos espacial

presencia <- st_as_sf(
  presencia,
  coords = c("decimalLongitude", "decimalLatitude"),
  remove = FALSE, # conservar las columnas de las coordenadas
  crs = 4326
)

# Gráfico ggplot2
grafico_ggplot2 <-
  presencia |>
  st_drop_geometry() |>
  group_by(year) |>
  summarize(n = n()) |>
  ggplot(aes(x = year, y = n)) +
  geom_line() +
  geom_point(
    aes(
      text = paste0(
        "Año: ", year, "\n",
        "Cantidad de registros: ", n
      )
    )
  ) +
  ggtitle("Cantidad de registros de presencia por año") +
  xlab("Año") +
  ylab("Cantidad de registros de presencia") +
  labs(caption = "Fuente: GBIF") +
  theme_economist()

Warning in geom_point(aes(text = paste0("Año: ", year, "\n", "Cantidad de
registros: ", : Ignoring unknown aesthetics: text

# Gráfico plotly
ggplotly(grafico_ggplot2, tooltip = "text") |> 
  config(locale = 'es')

# Gráfico ggplot2
grafico_ggplot2 <-
  presencia |>
  st_drop_geometry() |>
  group_by(year) |>
  summarize(n = n()) |>
  ggplot(aes(x = year, y = n)) +
  geom_line() +
  geom_point(
    aes(
      text = paste0(
        "Año: ", year, "\n",
        "Cantidad de registros: ", n
      )
    )
  ) +
  ggtitle("Cantidad de registros de presencia por año") +
  xlab("Año") +
  ylab("Cantidad de registros de presencia") +
  labs(caption = "Fuente: GBIF") +
  theme_economist()

Warning in geom_point(aes(text = paste0("Año: ", year, "\n", "Cantidad de
registros: ", : Ignoring unknown aesthetics: text

# Gráfico plotly
ggplotly(grafico_ggplot2, tooltip = "text") |> 
  config(locale = 'es')

# Mapa
leaflet() |>
  addTiles(group = "Mapa general") |>
  addProviderTiles(
    providers$Esri.WorldImagery, 
    group = "Imágenes satelitales"
  ) |>
  addProviderTiles(
    providers$CartoDB.Positron, 
    group = "Mapa blanco"
  ) |>  
  addCircleMarkers(
    # capa de registros de presencia (puntos)
    data = presencia,
    stroke = F,
    radius = 3,
    fillColor = 'red',
    fillOpacity = 1,
    popup = paste(
      paste0("<strong>País: </strong>", presencia$country),
      paste0("<strong>Localidad: </strong>", presencia$locality),
      paste0("<strong>Fecha: </strong>", presencia$eventDate),
      paste0("<strong>Fuente: </strong>", presencia$institutionCode),
      paste0("<a href='", presencia$occurrenceID, "'>Más información</a>"),
      sep = '<br/>'
    ),
    group = "Registros de Bradypus variegatus"
  ) |>
  addLayersControl(
    baseGroups = c("Mapa general", "Imágenes satelitales", "Mapa blanco"),
    overlayGroups = c("Registros de Bradypus variegatus"))

# Consulta a WorldClim
clima <- worldclim_global(var = 'bio', res = 10, path = tempdir())

# Nombres de las variables climáticas
names(clima)

 [1] "wc2.1_10m_bio_1"  "wc2.1_10m_bio_2"  "wc2.1_10m_bio_3"  "wc2.1_10m_bio_4" 
 [5] "wc2.1_10m_bio_5"  "wc2.1_10m_bio_6"  "wc2.1_10m_bio_7"  "wc2.1_10m_bio_8" 
 [9] "wc2.1_10m_bio_9"  "wc2.1_10m_bio_10" "wc2.1_10m_bio_11" "wc2.1_10m_bio_12"
[13] "wc2.1_10m_bio_13" "wc2.1_10m_bio_14" "wc2.1_10m_bio_15" "wc2.1_10m_bio_16"
[17] "wc2.1_10m_bio_17" "wc2.1_10m_bio_18" "wc2.1_10m_bio_19"

# Definir la extensión del área de estudio
area_estudio <- ext(
  min(presencia$decimalLongitude) - 5, 
  max(presencia$decimalLongitude) + 5,
  min(presencia$decimalLatitude) - 5, 
  max(presencia$decimalLatitude) + 5
)

# Recortar las variables bioclimáticas al área de estudio
clima <- crop(clima, area_estudio)

# Paleta de colores de temperatura
colores_temperatura <- colorNumeric(
  # palette = "inferno",
  # palette = "magma",
  palette = rev(brewer.pal(11, "RdYlBu")),
  values(clima$wc2.1_10m_bio_1),
  na.color = "transparent"
)

# Paleta de colores de precipitación
colores_precipitacion <- colorNumeric(
  # palette = "viridis",
  # palette = "YlGnBu",  
  palette = "Blues",
  values(clima$wc2.1_10m_bio_12),
  na.color = "transparent"
)

# Mapa
leaflet() |>
  addTiles(group = "Mapa general") |>
  addProviderTiles(
    providers$Esri.WorldImagery, 
    group = "Imágenes satelitales"
  ) |>  
  addProviderTiles(
    providers$CartoDB.Positron, 
    group = "Mapa blanco"
  ) |>
  addRasterImage( # capa raster de temperatura
    clima$wc2.1_10m_bio_1,
    colors = colores_temperatura, # paleta de colores
    opacity = 0.6,
    group = "Temperatura",
  ) |>
  addRasterImage( # capa raster de precipitación
    clima$wc2.1_10m_bio_12,
    colors = colores_precipitacion, # paleta de colores
    opacity = 0.6,
    group = "Precipitación",
  ) |>
  addCircleMarkers(
    # capa de registros de presencia (puntos)
    data = presencia,
    stroke = F,
    radius = 3,
    fillColor = 'red',
    fillOpacity = 1,
    popup = paste(
      paste0("<strong>País: </strong>", presencia$country),
      paste0("<strong>Localidad: </strong>", presencia$locality),
      paste0("<strong>Fecha: </strong>", presencia$eventDate),
      paste0("<strong>Fuente: </strong>", presencia$institutionCode),
      paste0("<a href='", presencia$occurrenceID, "'>Más información</a>"),
      sep = '<br/>'
    ),
    group = "Registros de Bradypus variegatus"
  ) |>  
  addLegend(
    title = "Temperatura",
    values = values(clima$wc2.1_10m_bio_1),
    pal = colores_temperatura,
    position = "bottomleft",
    group = "Temperatura"
  ) |>
  addLegend(
    title = "Precipitación",
    values = values(clima$wc2.1_10m_bio_12),
    pal = colores_precipitacion,
    position = "bottomleft",
    group = "Precipitación"
  ) |>  
  addLayersControl(
    # control de capas
    baseGroups = c("Mapa general", "Imágenes satelitales", "Mapa blanco"),
    overlayGroups = c("Temperatura", "Precipitación", "Registros de Bradypus variegatus")
  ) |>
  hideGroup("Precipitación")

coordenadas_presencia <- data.frame(
  decimalLongitude = presencia$decimalLongitude,
  decimalLatitude = presencia$decimalLatitude
)

coordenadas_presencia <- unique(coordenadas_presencia)

# Dividir los datos en entrenamiento y evaluación
set.seed(123)
n_presencia <- nrow(coordenadas_presencia)

indices_entrenamiento <- sample(
  1:n_presencia, 
  size = round(0.7 * n_presencia)
)

entrenamiento <- coordenadas_presencia[indices_entrenamiento, ]
evaluacion <- coordenadas_presencia[-indices_entrenamiento, ]

# Los datos de clima deben convertirse al formato que usa el paquete raster
# debido a es este el que acepta el paquete dismo
clima_raster <- raster::stack(clima)

# Ejecutar el modelo
#modelo_maxent <- maxent(x = clima_raster, p = entrenamiento)
modelo_maxent <- domain(x = clima_raster, p = entrenamiento)

# Crear un raster basado en el modelo
prediccion <- predict(modelo_maxent, clima_raster)

# Extraer valores del raster de predicción para los puntos de evaluación de presencias
eval_pres <- terra::extract(
  prediccion, 
  evaluacion[, c('decimalLongitude', 'decimalLatitude')]
)

# Generar puntos aleatorios de ausencia
ausencias <- randomPoints(mask = clima_raster, n = 1000)

# Extraer valores del raster de predicción para los puntos de ausencias
eval_aus <- terra::extract(
  prediccion, 
  ausencias
)

# Generar estadísticas de evaluación del modelo
evaluation <- evaluate(p = eval_pres, a = eval_aus)

# Graficación de la curva ROC
plot(evaluation, 'ROC')

# Paleta de colores del modelo
colores_modelo <- colorNumeric(
  palette = c("white", "black"),
  values(prediccion),
  na.color = "transparent"
)

# Mapa
leaflet() |>
  addTiles(group = "Mapa general") |>
  addProviderTiles(
    providers$Esri.WorldImagery, 
    group = "Imágenes satelitales"
  ) |>  
  addProviderTiles(
    providers$CartoDB.Positron, 
    group = "Mapa blanco"
  ) |>
  addRasterImage( # capa raster de temperatura
    clima$wc2.1_10m_bio_1,
    colors = colores_temperatura, # paleta de colores
    opacity = 0.6,
    group = "Temperatura",
  ) |>
  addRasterImage( # capa raster de precipitación
    clima$wc2.1_10m_bio_12,
    colors = colores_precipitacion, # paleta de colores
    opacity = 0.6,
    group = "Precipitación",
  ) |>
  addRasterImage( # capa raster del modelo de distribución
    prediccion,
    colors = colores_modelo,
    opacity = 0.6,
    group = "Modelo de distribución",
  ) |>  
  addCircleMarkers(
    # capa de registros de presencia (puntos)
    data = presencia,
    stroke = F,
    radius = 3,
    fillColor = 'red',
    fillOpacity = 1,
    popup = paste(
      paste0("<strong>País: </strong>", presencia$country),
      paste0("<strong>Localidad: </strong>", presencia$locality),
      paste0("<strong>Fecha: </strong>", presencia$eventDate),
      paste0("<strong>Fuente: </strong>", presencia$institutionCode),
      paste0("<a href='", presencia$occurrenceID, "'>Más información</a>"),
      sep = '<br/>'
    ),
    group = "Registros de Bradypus variegatus"
  ) |>  
  addLegend(
    title = "Temperatura",
    values = values(clima$wc2.1_10m_bio_1),
    pal = colores_temperatura,
    position = "bottomleft",
    group = "Temperatura"
  ) |>
  addLegend(
    title = "Precipitación",
    values = values(clima$wc2.1_10m_bio_12),
    pal = colores_precipitacion,
    position = "bottomleft",
    group = "Precipitación"
  ) |>
  addLegend(
    title = "Modelo de distribución",
    values = values(prediccion),
    pal = colores_modelo,
    position = "bottomright",
    group = "Modelo de distribución"
  ) |>  
  addLayersControl(
    # control de capas
    baseGroups = c("Mapa general", "Imágenes satelitales", "Mapa blanco"),
    overlayGroups = c(
      "Temperatura",
      "Precipitación",
      "Modelo de distribución",
      "Registros de Bradypus variegatus"
    )
  ) |>
  hideGroup("Temperatura") |>
  hideGroup("Precipitación")

Warning in colors(.): Some values were outside the color scale and will be
treated as NA